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FOREWORD 

The  estimation  of  the  time  delay  between  signals  arriving  at  two  spatially 
separated  sensors  is  investigated,  for  the  purpose  of  target  localization. 

The  real  data  were  obtained  during  a  sea  test. 


F.  B.  SANCHEZ  ~0 
By  direction 
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INTRODUCTION 


This  report  examines  at  first  the  estimation  of  time  delay  r (t)  and  secondly 
the  estimation  of  time  delay  and  delay  rate  r  (t)  (doppler  shift),  between 
signals  received  at  two  spatially  separated  sensors  in  the  presence  of 
uncorrelated  noise.  The  mathematical  models  for  both  estimation  procedures 
assume  stationarity  for  the  characteristics  of  signal  and  noise  for  a  finite 
observation  time  T.  It  is  also  assumed  that  the  spectral  characteristics  of  the 
signal  are  known. 

ONE  DIMENSIONAL  SEARCH  FOR  THE  ESTIMATION  OF  THE  TIME  DELAY  r (t) 

FORMULATION  AND  ANALYSIS. 


A  signal  s(t)  emanating  from  an  acoustic  source  and  received,  in  the 
presence  of  noise,  at  two  spatially  separated  sensors  can  be  mathematically 
modeled,  as  shown  in  Figure  1,  as  follows: 


xi(t)  -  s(t)  ♦  ni(t) 


.(1-1) 


x2(t)  -  s(t-r(t))  +  n2(t) 


FIGURE  1  SIGNAL  DELAYED  AND  CORRUPTED  BY  NOISE, 
RECEIVED  AT  THE  TWO  SENSORS 


Xt(t) 


Xjjlt) 


(1-2) 


It  is  assumed  that  s(t),  ni(t),  n2(t)  are  jointly  stationary  random 
processes,  and  the  signal  s(t)  is  uncorrelated  with  noise  n^(t)  and  n2(t), 
1*6*  ) 

E{s(t)ni*(t)f  -  E|s(t)n2*(t)}  -  0, 

where  E  denotes  expectation  and  *  complex  conjucate.  We  also  assume  that  we 
have  real  random  processes  and  furthermore, 
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E js(t)s*(t) }  =  Q,  E|ni(t)iij(t)|  =  Rfoj.i,  j  *  1,2 

The  time  delay  r  is  the  parameter  to  be  estimated  by  r'.  One  method  of  such 
estimation  (  [1]  )  is  to  maximize  the  crosscorrelation  function. 

Rxl’  x2^  *  EIx1^c^x2*  (1-3) 

The  argument  r ■'r  that  achieves  the  above  maximum  is  the  desired  time  delay 
estimate. 

For  the  given  model  we  have,  under  the  assumptions: 

x2<'>  -  Rs,s<r-^>  +  n2<^  <!“*> 


The  Fourier  transform  of  (1-4)  gives  the  cross-spectrum: 
Sx^,x2^  “  Ss,s^e  2irj^r+Sn^>n^(f ) 


Since 


Rn^n^r)  3  0.  then  ,  ^(f)  *  0 


From  (1-6),  equation  (1-5)  is  equivalent,  in  the  time  domain,  to: 
Rx  v  (t)  ■  R  (t)»5(t-t) 


X1,X2 


s  ,s 


(1-5) 


(1-6) 


(1-7) 


where  •  denotes  convolution. 

Equation  (1-7)  says  that  the  delta  function  has  been  "smeared"  by  the  signal. 

If  the  signal  s(t)  is  white  noise,  only  then  the  delta  function  is  not  spread 
and  the  r  is  the  position  of  the  delta  function.  For  this  special  case,  see 
Figure  7.  In  this  figure,  the  white  noise  signal  received  at  one  sensor 
relative  to  the  other,  is  delayed  by 

T*  0,  "1/8,  “1/4  seconds. 

Each  correlogram  is  averaged  for  4  seconds.  During  the  first  24  seconds  we  have 
infinite  signal  to  noise  ratio,  while  for  the  remaining  time  we  have  S/N  ratio  * 
-15  dB.  Different  delay  rates,  positive  and  negative,  are  also  illustrated. 

IMPLEMENTATION. 

The  above  time  delay  estimation  procedure  was  implemented  in  the  MAP-300 
array  processor  and  the  POP  11/34  as  the  host  computer.  As  a  continuation  of 
the  block  diagram  in  Figure  1,  consider  the  block  diagram  of  Figure  2, 
describing  the  algorithm. 
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f„  SAMPLING  FREQUENCY 


*1 


*2 

Xl 


x2 


FIGURE  2  BLOCK  DIAGRAM  OF  THE  TIME  DELAY  ESTIMATOR 


FIGURE  3  BLOCK  DIAGRAM  OF  THE  TIME  DELAY  ESTIMATOR 
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From  Figure  2,  the  time  sequences  x^(t),  X2(t)  go  through  the  low-pass 
filters  H^(f),  ^(f).  The  maximum  signal  frequency,  fmax,  allowed  through 
these  filters  is  such  that  fs  ^  2.5fmax>  where  fs  is  the  sampling 
frequency.  Once  this  is  done,  the  filtered  sequences  xj^(t),  X2(t)  enter  the 
ADAM  (Analog  Digital  Acquisition  Module)  where  they  are  sampled  at  the  rate  fs, 
and  multiplexed.  This  process  is  controlled  by  the  software.  Once  sampled,  one 
of  the  double  buffers  labeled  T1MEA,  and  containing  both  digitized  sequences 
xi(t),  X2(t),  is  processed  by  the  MAP-300,  while  the  other  double  buffer 
TIMEB  is  being  filled  with  the  new  digitized  data.  This  guarantees  that  no  data 
is  lost  during  the  processing  and  furthermore,  it  guarantees  parallel  processing 
between  the  ADAM,  the  MAP-300  and  the  host  computer. 

Figure  3  shows  the  demultiplexing  of  the  time  sequences,  which  are  Chen 
FFT'd,  cross-multiplied  after  one  has  been  conjucated  and  the  cross  product 
inverse  FFT'd.  From  the  resulting  complex  time  sequence,  we  obtain  a  real  time 
sequence  by  taking  the  power  in  each  of  the  cells  of  this  complex  sequence.  The 
resulting  real  time  sequence  is  the  correlogram  (Eq.  (1-7)),  which  is  averaged 
for  T*  seconds.  The  averaged  correlograms  are  then  plotted  on  the  Tektronix 
4014  by  using  Fortran  calls  from  the  host  computer. 

It  is  of  importance  to  note  that  all  of  the  above  procedure  depicted  in 
Figure  3  is  done  in  the  MAP  300  in  such  a  way  as  to  achieve  maximum  speed  during 
the  real  time  processing  of  the  data.  It  is  also  important  to  note  the 
sychronization  (parallel  processing)  between  the  MAP-300  and  the  host  computer. 
The  time  during  which  it  takes  Che  host  to  plot  the  (k-l)ch  correlogram,  the 
MAP-300  is  processing  che  data  obtained  from  the  double  buffer  TIMEB  (part  (B) 
of  Figure  2).  So  by  che  time  the  plotting  of  the  (k-l)ch  correlogram  is 
completed,  the  kth  correlogram  is  available  to  the  host  for  plotting. 

The  duration  T,  of  the  FFT  shown  in  Figure  3  is  T  *  N/ fs ,  where  N  is  the 
number  of  points  in  each  of  the  channels  xj^lt),  X2(t),  and  fs  is  the 
sampling  frequency.  A  typical  sampling  duration  T  is  T  s  1  sec.  This  provides 
the  capability  to  search  for  r"  in  a  delay  window  from  -1/2  to  +1/2  seconds. 

Based  on  the  distance  between  the  sensors  the  delay  window  can  vary  by 
appropriately  varying  N  or  fs  or  both.  Also,  the  number  of  seconds  for  the 
post-integration  of  the  correlogram  is  a  variable.  If  T'  *  4  sec  and  T  *  1  sec 
(duration  of  a  time  cut)  and  if  we  have  10  minutes  of  data  to  process  (600  time 
cuts),  the  number  of  correlograms  plotted  is  N^=  600/T'  *  150. 

Following  a  spectral  analysis  of  the  signal,  based  on  its  spectral 
characteristics,  we  may  want  to  process  only  a  spectral  window  a round  some 
center  frequency  f0.  The  outlined  algorithm  is  capable  of  doing  this,  by 
extracting  any  spectral  window  for  every  channel  and  searching  for  r,  based  on 
the  energy  of  the  signal  inside  this  spectral  band. 

The  algorithm  described  above  is  successful  in  estimating  the  time  delay 
between  signals  received  at  sensors  closely  spaced  together  and  between  signals 
from  a  nonmaneuvering  target  moving  at  a  relatively  low  speed. 


IfCnapp  C.  H.  and  Carter  G.  C. ,  "The  Generalized  Correlation  Method  for  Estimation 
of  Time  Delay,"  IEEE  Trans.  Acous.,  Speech,  and  Signal  Proc.,  Vol.  ASSP-24,  No. 
pp.  320-327,  Aug  1976. 
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It  is  known  that  the  received  frequency  at  a  sensor  due  to  the  relative 
motion  between  the  target  and  sensor  is  given  by 


f  3  f  <1  -  XS2S1  >  U-8) 

c 

where: 

f:  received  frequency 
f0:  radiated  frequency 
v:  target  velocity 
c:  speed  of  underwater  sound 

0:  angle  between  target  velocity  vector  and  the  line  between 
sensor  and  target. 

Equation  (1-8)  shows  that  it  is  necessary  to  compensate  for  the  high 
doppler  shift  arising  from  maneuvering  targets  or  targets  moving  at  high 
speed.  The  analysis  and  algorithm  that  follow  consist  of  a  two  dimensional 
search  in  which  in  addition  to  the  estimation  of  time  delay  r (t)  an  estimate  of 
the  time  delay  rate  f(t)  is  obtained  to  accommodate  for  the  cases  mentioned 
above . 


TWO  DIMENSIONAL  SEARCH  FOR  THE  ESTIMATION  OF  THE  TIME  DELAY  r(t)  AND  TIME  DELAY 
RATE  f  (t) 


FORMULATION  AND  ANALYSIS. 

As  in  Part  1  of  this  report,  let  s(t)  be  the  radiated  signal  and  let 
n^(t),  n2(t)  the  noise  received  at  the  two  sensors.  Under  the  same 
assumptions  we  have 

x].(t)  *  s(t)  +  n^(t)  (2-1) 

X2(t)  -  s(t-r+  at)  +  n2  (t)  (2-2) 

where  r  and  a  are  the  delay  and  delay  rate  respectively  to  be  estimated. 
Furthermore,  let  the  time  series  x^(t),  X2(t)  be  heterodyned  and 
appropriately  filtered  with  weights  w(t).  The  block  representation  of  this 
procedure  is  shown  in  Figure  4. 

««>  -  /CN  t{t) 

T) 

•  <t) 

FIGURE  4  HETERODYNE  AND  LOW-PASS  OF  THE  TIME  SEQUENCE 
We  have  z(c)  a  x(t)  e(t) 

and  y(t)  ■  z(t)^(t)  *  x(t)e(t)®v(t)  (2-3) 


w<t) 


V(t> 
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where  e(t)  “  e”iwhc 


w(t)  m  (1-cos^-),  0  <  t  <  T 


Che  Hanning  weights. 

For  each  time  series  xi(t),  X2(t)  we  have: 

00 

yi(t)  */«,xi(si  +  t)e“^,h^sl  +  c^w(si)dsi 
00 

y2(t)  *  /  X2(s2  +  t)e_i‘*'h(s2  +  c)w(s2)ds2 
00 

But  X2(t)  =*  xi(t-r0  -K»t)  ■  xi((l+ar)t-r0) 

Also  let  E  |xi(c)x*  (C’)|  *  R(c-t') 

Then  Che  crosscorrelation  function  of  the  output  series 
yi(t),  y2(t)  is 

R  (t)  =■  E|y.(t-T)y2  (t)  \  m 

yuyz  m  L 


(2-4) 


(2-5) 

(2-6) 


(2-7) 


(2-8) 


r°°  m 

I-aa  Xao  X1  ( s  1+ 1~  r ) X2  ( s 2+t )  | e~ ^h  ^ s  1  _s2 ~ T (  s  i  )w(s2)ds^ds2 

CO  oo 

■lot  —oo  E{*i(*l+t-  r)x*  (  (1+Cl)  (s2  +  t  )-  t  q)  \  e~4*,h^sl-s2'~1’  )w(si)w(s2)ds^ds2 

00  CO 

J  J*  R[  (s^  +  t-r  )-(  ( 1+a)  (s2+t )-  r0)  ]  e'^h  ^sl_s2~r  ^w(s  i  )w(s2  )dsids2 


If  the  signal  is  white  noise  then 


Rx^>x^  *  E|  xi(t)x*(t')}  *  R(t-t')  *  5(t-t') 
and  Eq.  (2-9)  becomes  for 

s  i  *  (1+a  )s2  +  (r-7o^  +ac  and  s2  a  s » 

CO 

Rv  „  (r)  *  /  e-*^  (<*s+at-  ^  )w(s  )w(  (l+o)s  +  (  r-  r  +at))ds 

j  l  'j  2  — 00 


(2-9) 


(2-10) 


(2-11) 

(2-12) 


We  can  write  r-r  +a t  as 


r  -  T0+ac  *  «]  +  e2t  where 
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* ^  is  Che  error  in  estimating  rQ(  the  time  delay) 

and  e2  is  Che  error  in  estimating  a (the  time  delay  rate) 

IMPLEMENTATION. 

The  implementation  of  equations  (2-5)  and  (2-6)  in  the  MAP-300  was  done  by 
using  the  DCVM  (Discrete  Convolution  Multiply)  and  CXMUL  (Complex  Multiply) 
calls.  If  we  let  T  be  the  sampling  time,  f^  the  heterodyne  frequency,  fs 
the  sampling  frequency  and  N  the  number  of  points  in  the  time  series,  then  the 
"select"  function  e(t)  in  the  discrete  case  becomes: 


e.  3  e^^ht 
k  k 


g-yrifhkT/N  =  e-2n.fh(T/N)k 


e-2*l(fs)k 


So,  from  equation  (2-5)  or  (2-6) 


(2-13) 


y 


k 


wk®xkek 


M-l  £h  M  M 

2  w(m)e~2?rifs  (®+k4)  x(m+k2T) 

m“0 


y 


k 


e-ikC 


Lh  Mt. 

f  2  ; 
s 


M-i  fh 

2  w(m)e-2  x(m+k~) 

m*0  H 


(2-14) 


M 

where  4  is  Che  decimation  factor,  corresponding  to  a  4:1  redundancy  in  the 
weighting  of  the  time  series  x^(t),  X2(t)  by  the  Hanning  weight  coefficients 
given  by 


w(m)  *  4  (1  -  cos-^) ,  m“0,l . M-l. 

2  M 

To  determine  the  frequency  response  of  the  heterodyned  and  low-passed  time 
series  x(t),  let 


w(t)  -  4  (1  +  cosili)  for  -T/2  <  t  <  T/2 
2  T 


let 


e(t)  *  e-2ir*-^ht 
and  let  x(t)  *  e2irift>  Then  from  Figure  4  we  have: 

y(t)  “4  /  (1  +  cos4p-  J  e~2 'rifh  ^c-s  ^x( t-s  )ds 

2-T/ 2  T 


(2-15) 
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y(t) 


7  Tf2  [1+4  (eiYs  +  e-iYs)]e27ri(f-fh)t  e-2iri(f-fh  >sds 
2  -T/ 2  2 


(2-16) 


where  T  »  .  If  5f  =  (f-fp,  then 


y(t)  ■  y-  e2,ri($f)t  Tf22e-2,ri^f  )s+e-2xi(5f  +  -)3+e-27ri(5f  -  -)sds 
4  _T/2  T  T 

In  discrete  form  (2-17)  becomes: 


(2-17) 


M.  1 
2 


„  -PP,  2  „  <5f  if  1  <5£  1 

4e-2^i(-^)k  ,2  [2e-27ri(— )m  +e-2n(f  +  pm+e-2»-i(—  ~  M>ml 


8  'k  m  M 

““‘I 


rZ2  ££ 

It  can  be  shown  that  J  e-2*i.(f^)s  ds  *  ^*0 


1 

-T/2 


sin(0rT) 


(2-18) 

(2-19) 


£f 

where  =*  (f  ).  Using  (2-19),  Eq.  (2-18)  can  be  approximated  by, 
s 


4e-2*'i0ky  _  2sin(7r0M)  sin7r(fi+  ft)M  sinTrO-  H)M 

k  "  *0  +  tt(/3+  p  +  "-(0-  pM 

M  M 


(2-20) 


or 


4e-2iri0ky  2s  in  Orff  M)  s  in (*0  M)  sinOr0M) 

k  *0  jr(0+  p  tt(0  —  p 


(2-21) 


Eqivalently: 


4e-2iri£ky  sinGrflM)  t2  _  1 _ L___) 

k  3  TT  0  (j3+  p  (0-  p 


(2-22) 


4jre~27ri^ky  m  _ M _  sin&r0M) 

k  0(02  -  p) 

9 

M" 

P 

2re“2’ri5ky  _  jV*2 _  sinGr0M) 

e<  i  - e2> 


(2-23) 


(2-24) 
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*■  *■ 

Letting  r  ■  M/3  *  M(— — )  in  (2-24)  we  get: 

s 

2ire“2,ri0k  y  „  jj  .  sin(irr) 
k  r(l-r2) 

or 

2e  2ir^k  y  _  sin(irr)  „  -sin[x(r-l)]  (2-25) 

M  k  *r(l-r2)  rr(l-r2) 

The  right-hand  ride  of  equation  (2-25)  is  the  frequency  response  Y(r)  of  the  time 
series 

8 

Clearly  Y(0)  -  1,  Y(.5)  -  ,  Y(l)  -  .5. 

3ir 


M(f-fh) 

So  for  r  *  1  ■  — - -  we  have, 

s 

at  6-dB  down,  Che  6-dB  bandwidth 


1  -  <f-*h> 


f 

s 

M 


or 

2f 

®6-dB  “  "/ 


(2-26) 


The  block  diagram  describing  the  implementation  of  the  two  dimensional  search 
algorithm  is  shown  in  Figure  5.  It  is  a  continuation  of  part  (A)  of  the  one 
shown  in  Figure  2,  in  the  first  part  of  the  report. 

As  shown  in  Figure  5,  the  double  buffer  containing  the  digitized  time  series 
xj(t)  and  X2(t)  is  demultiplexed  by  the  software  into  two  channels  xi(t), 

X2(t)  representing  Che  received  signal  s(t)  at  the  two  sensors,  of  respective 
frequencies  fQ  and  f0(l+x),  where  t  is  the  doppler  shift  to  be 
estimated.  The  two  time  sequences  xi(t),  X2(t)  are  then  heterodyned.  The 
effect  of  the  heterodyning  is  to  select  the  spectral  window  of  interest  around 
some  center  frequency  fj, .  The  width  of  this  spectral  window  is  approximately 
given  by  equation  (2-26).  As  a  result  of  the  low-pass  filtering  and 
heterodyning,  the  output  complex  time  sequences  yi(t),  y2(t)  have  respective 
frequencies  (fQ  -ffo)  and  (f0  +  f0r  -  f^).  At  this  point  yi(t), 
y2(t)  are  the  input  channels  to  the  crosscorrelator  that  will  estimate  r(t) 
and  f ( t ) . 


Let  N  *  number  of  complex  points  to  be  processed  by  the  crosscorrelator. 
(*size  of  yi(t),  y2(t)). 

Let  M  *  number  of  Hanning  weight  coefficients. 
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r 


CROSSCORRELATOR 


f0r  (beat  frequency) 
rttime  delay) 


FIGURE  5  BLOCK  DIAGRAM  OF  THE  TIME  DELAY  AND  TIME  DELAY  RATE 
ESTIMATOR 
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Because  of  Che  4:1  redundancy  (M/4  decimacion  factor  in  Che  convolution, 

Eq.  (2-14))  in  the  weighting  of  the  time  sines  x^(t),  X2(t),  we  ®ust  have 
MN/4  points  for  each  xi(t),  X2(t).  If  fs  is  the  sampling  frequency,  then 

MN 

4 

At  * 

s 

gives  the  number  of  seconds  corresponding  to  the  complex  time  series  y^lt), 
y2(t)  entering  the  crosscorrelator.  The  number  of  points  FFT'd  is 
N'fft  *  2N  (we  have  a  2:1  redundancy  in  Che  FFT).  Typical  values  for  the 
parameters  defined  above  are:  M  *  128  pts,  N  *  64  pts,  fs  ■  1024  pts/sec. 
These  values  for  equation  (2-26),  give 


2f 
_ s 

II 


16Hz 


frequency  window  to  be  crosscorrelated.  Furthermore  the  duration  At  of  the  time 
cuts  (size  of  yi(t),  y2(t))  entering  the  crosscorrelator  is, 


MN 


At 


4  (32)  (64) 

"2 —  *  ~ —  ■  2  seconds. 


1024 


Moreover,  the  FFT  sampling  frequency,  (decimated  frequency),  fsFFT*  *s 


fSFFT*  ■  -y-  *  32pts/sec, 


per  channel.  Therefore  the  FFT  time 
N'FFT  2N 


.  tFFt  -  f  -  N 

SFFT  -n- 
4t 


2At  ■  4  seconds 


because  of  the  2:1  redundancy  in  the  FFT.  The  FFT  resolution  is 


rFFT 


FFT 


2At 


Hz/bin. 


Figure  6  shows  how  the  2:1  redundancy  in  the  FFT  and  the  variable  delay 
scheme  are  implemented  in  the  MAP-300  array  processor.  For  this  algorithm  we 
have  chosen  a  slow  time  delay  search  (through  shifting  in  the  time  domain)  and  a 
fast  time  delay  rate  search  (via  an  FFT). 


15 


NSWC  TR  81-303 


W(t):  HANNING  WINDOW 


CHANNEL  1. 


r 

* 


A1 


A2 


1 


A3 


* 


\t*1 

\ 

\ 

\ 


t-2 


I  / 

CHANNEL  2.T  I  • 


t-0  'Uf. 


B1 


B2 

»■  - 


t-1 

SHIFT 


t-2 


>* 

I  / 

I  / 

^  B3 

I 

I 

I - 


/  \t-3 

/  \ 

% 


B4 


\  I 
\  I 


t-3 

SHIFT 


t-4 


t-4 


NEW  DATA 
yi<t) 

•tt-5 


NEW  DATA 

y2<*> 

att-5 


FIGURE  6  VARIABLE  TIME  DELAY  SCHEME  AND  502  OVERLAP 
IN  THE  WEIGHTING  OF  CHANNEL  1 


Segments  A1  through  A3  and  B1  through  B4,  each  contain  the  N  points  entering 
the  crosscorrelator.  The  variable  t  represents  units  of  time.  In  this  case 
each  unit  of  time  is  of  duration  At  seconds.  For  the  values  considered 
previously,  N-64  pts  and  At  -  2  seconds.  Before  any  new  data  is  received, 
segments  A1  and  A2  are  weighted  with  a  2N-point  Hanning  weight.  The  resulting 
vector  of  size  2N  is  V^.  For  channel  2,  to  accomodate  for  a  delay  window  of 
♦  r(t)  seconds,  we  pick  the  vector  V2 ,  which  starts  ^  points  to  the  left  of 
t-1,  where 


N  r(c) 

\  “  Ae  P's- 

We  then  form  the  point-by-point  vector  product  U  -  V*  V2,  which  is 
FFT’d.  Following  this  we  form  a  new  vector  V2  of  length  2N  points,  starting 
at  (r^  -1)  points  to  the  left  of  t  ■  1  and  again  we  form  a  new  U  vector  which 
is  again  FFT’d.  This  procedure  is  repeated  (2t[<  +  1)  times  so  that  we  sweep 
through  all  points  representing  the  negative  and  positive  delay  window.  Once 
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Che  delay  search  is  completed  for  this  fixed  time  cut  of  duration  At  seconds,  we 
shift  Che  segments  A2-»A1,  A3-A2 ,  B2-»B1,  B3-*B2 ,  B4-*B3  and  the  newest  segment  of 
data  for  y^d),  y2^t)  enters  in  A3  and  B4  respectively,  in  such  a  way  that 
the  data  points  are  ordered  in  time.  Moreover,  Che  units  of  time  are  also 
updated  by  1.  The  delay  search  is  then  repeated  for  Che  new  time  cut.  This 
procedure  enables  us  to  search  for  r (t)  within  a  +At  *  +2  second  window. 

So  for  every  time  cut  and  for  every  delay  point  r^,  after  y^(t)  has  been 
Hanning  weighted  in  a  2:1  fashion,  conjucated  and  multiplied  by  a  delayed 
y2(t),  we  take  the  FFT  and  calculate  the  power  spectrum  of  the  product.  We 
then  search  for  the  location  of  the  peak  in  the  power  spectrum,  from  which  we 
can  get  an  estimate  of  the  delay  rate  r  (see  below).  With  an  appropriate 
shifting  of  the  power  spectrum  the  zero  delay  rate  (f«  0)  corresponds  to  the 
frequency  bin  in  the  middle  of  the  power  spectrum.  This  enables  us  to  look  only 
at  a  window  centered  at  the  middle  of  the  power  spectrum.  The  width  of  the 
window  can  vary  as  a  function  of  Che  f  chat  is  to  be  estimated. 

The  change  Ar  in  Che  delay  T  ,  per  update  of  duration  At,  is  given  by 


At  -  (At)r 


2  tSFFT* 


1  2N  •  N  • 

2  f  T“  f  T 

SFFT  S FFT 


where  *  is  the  delay  rate.  The  frequency  cell  where  the  peak  in  the  power 
spectrum  shows  up  for  every  time  cut  (update)  and  for  every  delay  point  T ^  is 
given  by, 


a  •  *  9V 

nFFT  *  for  TFFT  *  fhr  TFFT  *  *hr  ? -  (2~27 

SFFT 

The  position  of  the  peak  at  the  (Tfc)ch  delay  point  is 

nTl  *  fSFFTT 
k 


So  for  fsFFT  *  32pts/sec  and  r*  1/16  sec, 

n  -  2 
rk 


positions  from  the  zero  delay  position. 

The  rate  An  at  which  the  position  of  the  peak  in  the  power  spectrum  is 
moving  for  every^pdate  At  is 


An 

rk 


|  (2N)  f  «  Nf 


(2-28 
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^  , 

Tk 


Nr 


f$FFT 

N - n 

2Nfh 


FFT 


1  ,  FSFFT 

2  ( — : — m)  nFFT 

fh 


where  is  Che  heterodyne  frequency. 
ILLUSTRATION  OF  RESULTS 


The  figures  Chac  follow  illustrate  some  results  obtained  using  the 
algorithms  described  in  this  report.  In  all  the  figures  the  locations  of  three 
or  four  strongest  peaks  in  the  correlogram  were  plotted  by  the  host  computer. 

Figure  8  shows  the  time  delay  track  of  a  signal  emanating  from  a  target  of 
interest.  This  track  was  obtained  using  the  one  dimensional  sheme  with  no  doppler 
compensation.  The  spacing  between  the  sensors  is  approximately  250  yards.  We 
estimated  that  the  target's  speed  is  about  6  knots  and  that  the  target  is 
located  at  about  2600  feet  from  one  of  the  sensors  at  the  CPA  time  to  this 
sensor.  Ws  processed,  for  20  minutes,  a  spectral  window  of  width  256  Hz. 

Tne  scenario  for  the  time  delay  tracks  shown  in  Figures  10  and  11  is  shown 
in  figure  9. 


11:11:46 


11:01:27 


FIGURE  9  SCENARIO  FOR  FIGURES  10  AND  11 

The  target  is  a  towed  sound  source  radiating  broadband  energy  from  200  to 
1000  Hz.  It  is  moving  at  a  speed  of  6  knots. 

Figure  10  shows  the  time  delay  Crack  between  a  Q57A  buoy  and  a  sensor  on  a 
Horizontal  Planar  Array.  The  distance  between  them  is  approximately  1500  feet. 
It  is  a  30  minute  run,  the  correlograms  are  averaged  and  plotted  every  4  seconds 
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KUV'IISCC 

Figure  10  Time  delay  track  between  a  Q57A  buoy  and  a  sensor  on  a  HPA 
Three  strongest  peaks  in  correloyrani 
Four  second  |jost  integrated  conelograms 
Cross  correlation  window  172  to  428  Hz 
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Flyura  11  Time  delay  back  between  two  sensors  on  a 
Ttiree  strongest  peaks  m  corretoyram 
Cross  correlation  window  200  to  800  Hi 
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and  the  three  largest  peaks  are  plotted.  For  the  first  10  minutes  a  direct 
signal  from  the  target  is  received  and  after  the  turn,  for  the  remaining  20 
minutes  the  bottom  bounce  signals  are  received. 

Figure  11  shows  the  time  delay  track  of  the  signal  from  the  target  of 
Figure  9,  received  at  two  sensors  located  on  a  Horizontal  Planar  Array. 


The  two  dimensional  search  algorithm  was  used  tor  the  scenario  depicted  in 
Figure  12. 
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FIGURE  12  SCENARIO  FOR  FIGURES  13  AND  14 


The  spacing  between  the  sensors  is  approximately  2200  yards  and  the  target 
is  moving  at  a  speed  of  about  15  knots. 

The  algorithm  was  successful  only  when  we  processed  frequency  windows  below 
40Hz.  Figure  13  and  14  show  the  time  delay  track  and  the  change  in  time  delay 
for  30  minutes.  For  both  cases  we  sampled  at  fs  *  512  pts/sec,  and  the 
decimated  frequency  was  16  Hz  (16  pts/sec  out  of  the  heterodyne  and  low-pass 
filter).  We  crosscorrelated  16  seconds  of  data  and  plotted  the  correlograms 
every  8  seconds.  The  time  delay  window  was  +2  seconds.  In  Figure  13  we 
processed  a  window  of  8Hz  centered  at  16Hz  and  in  Figure  14  the  same  frequency 
window  was  centered  at  30Hz.  As  we  increased  the  heterodyne  frequency,  the  time 
delay  track  dissappeared. 
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Figure  T3  Tim#  delay  track  with  do ppler  shift  compensation 
Heterodyne  frequency;  16  Hz 
Decimated  frequency;  16  Hz 
Correfogramt  every  8  seconds 
18  second  of  data  were  cross  correlated 


Figure  14  Tim*  dtlay  track  with  doppler  shift  compensation 
Heterodyne  frequency:  30  Hr 
Decimated  frequency  16  Hr 
Correlogram  every  3  seconds 
16  seconds  of  data  were  cross  correlated 
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SUMMARY  OF  RESULTS 


For  closely  spaced  sensors  and  slow,  non-maneuvering  targets  the  algorithm 
described  in  the  first  part  of  the  report  can  be  used.  In  all  other  cases,  when 
doppler  compensation  is  necessary,  the  two  dimensional  search  routine  can  be 
used,  although  its  success  is  limited  to  processing  broadband  energy  at  the  very 
low  end  of  the  spectrum. 
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